* FigureA5A6A7.do  12/07/23
* Alpert Powell Pacula 2018 semi-replication and then adding differences by age
* Figure A5.  Comparing adult and child opioid fatalities
* Sensitivity analysis for Figure A5

* (1) assemble dataset
cd ..
do Figure3prep
cd youthappendix
ren state CensusDivisionCDC
sort CensusDivisionCDC year
merge CensusDivisionCDC year using PolicyPanel
drop if year>2020
tab _merge
drop _merge
drop pop20042009 pop2010       // These are for the entire Census division, not just the age X gender therein
ren CensusDivisionCDC state

* (2) create regression variables
scalar youthfactor = 25
scalar baseyear = 2010
gen byte female = gender==`"{"Female"}"'
gen byte youth = age=="0-17"
encode state, gen(cdivnum)
gen deathrt = deaths*100000/pop
gen deathrtscaled = deathrt
replace deathrtscaled = deathrt*youthfactor if youth>0
gen lndeathrt = log(deathrt)
* (2b) declare national averages as scalars (for later use in lincom)
capture drop tmp
gen byte tmp = 0
preserve
tempfile navs
tempfile navmort
drop if age=="All"
keep if ucd=="Drug" & mcd=="Opioid"
gen deathrtscaledadult = deathrt if age=="18+"
gen deathrtscaledchild = deathrtscaled if age=="0-17"
gen lndeathrtadult = lndeathrt if age=="18+"
gen lndeathrtchild = lndeathrt if age=="0-17"
gen lndeathrtmen = lndeathrt if age=="18+" & female==0
gen lndeathrtboys = lndeathrt if age=="0-17" & female==0
collapse (median) tmp (mean) navdeathrtscaledchild=deathrtscaledchild navdeathrtscaledadult=deathrtscaledadult navlndeathrtchild=lndeathrtchild navlndeathrtadult=lndeathrtadult navlndeathrtboys=lndeathrtboys navlndeathrtmen=lndeathrtmen mean_initial_oxy=initial_oxy mean_havePDMPmust=havePDMPmust if age!="All" [aw=pop], by(year)
sort year
save `navmort', replace
tsset year
gen biann_initial_oxy = (mean_initial_oxy+l1.mean_initial_oxy)/2      // for use when aggregating adjacent years
gen biann_havePDMPmust = (mean_havePDMPmust+l1.mean_havePDMPmust)/2
reshape clear
reshape groups year 1999-2020
reshape cons tmp
reshape vars mean_initial_oxy mean_havePDMPmust biann_initial_oxy biann_havePDMPmust
reshape wide
sort tmp
save `navs', replace
restore
sort tmp
merge tmp using `navs'
tab _merge
drop _merge tmp
foreach yr of numlist 1999/2020 {
    scalar usapdmp`yr' = mean_havePDMPmust`yr'
	scalar usaoxy`yr' = mean_initial_oxy`yr'
    scalar bianusapdmp`yr' = biann_havePDMPmust`yr'
	scalar bianusaoxy`yr' = biann_initial_oxy`yr'
	drop mean_initial_oxy`yr' mean_havePDMPmust`yr' biann_initial_oxy`yr' biann_havePDMPmust`yr' // create scalars and drop the variables
}
sort year
merge year using `navmort'
tab _merge
gen navadultchildgap = navdeathrtscaledadult-navdeathrtscaledchild
gen navlnadultchildgapnog = navlndeathrtmen-navlndeathrtboys
drop _merge mean_initial_oxy mean_havePDMPmust navdeathrtscaledchild navlndeathrtadult navlndeathrtchild navlndeathrtmen navlndeathrtboys
sort cdivnum year
* (2c) APP 2018 "treatment variables" plus their interactions with youth
gen pretrend=(year-baseyear)*initial_oxy
gen posttrend=(year-2011)*initial_oxy*(year>2010)
gen shift=(year>2010)*initial_oxy
gen youthpre = youth*pretrend
gen youthpost = youth*posttrend
gen youthshift = youth*shift
gen youthoxy = youth*initial_oxy
gen PDMPmusty = youth*havePDMPmust
gen byte noreghavePDMPmust = 0
gen byte noregPDMPmusty = 0
gen byte noregposttrend = 0
gen byte noregshift = 0
gen byte noregyouthpost = 0
gen byte noregyouthshift = 0

local outstats noparen noaster adjr2 nobs addstat("s.e.",e(rmse))
local decs bdec(6) sdec(6) se rdec(4)

* (3) regression results: Reformulation
* (3a) closed spec to APP 2018, Table 3
reg deathrt i.cdivnum i.year female pretrend shift posttrend [aw=pop] if year>2007 & year<=2013 & age=="All" & mcd=="T400T401T404", cluster(cdivnum)
outreg2 using "FigureA5A6A7.outreg.xls", `decs' `outstats' addtext("Substance","Im opioids","Adults and minors","Pooled","Years","2008-2013","Groups excluded","None") excel replace
lincom shift+posttrend*2
* (3b) Kim sample ends 2016
reg deathrt i.cdivnum i.year female pretrend shift posttrend [aw=pop] if year>2007 & year<=2016 & age=="All" & mcd=="T400T401T404", cluster(cdivnum)
lincom shift+posttrend*2
* (3b) Kim uses earlier years too
reg deathrt i.cdivnum i.year female pretrend shift posttrend [aw=pop] if year<=2016 & age=="All" & mcd=="T400T401T404", cluster(cdivnum)
outreg2 using "FigureA5A6A7.outreg.xls", `decs' `outstats' addtext("Substance","Im opioids","Adults and minors","Pooled","Years","1999-2016","Groups excluded","None") excel append
lincom shift+posttrend*2
* (3d) just youth
reg deathrt i.cdivnum i.year female pretrend shift posttrend [aw=pop] if year>2007 & year<=2016 & age!="All" & youth>0 & mcd=="T400T401T404", cluster(cdivnum)
lincom shift+posttrend*2
* (3e) switching mortality indicator to all opioid (from just Im)
reg deathrt i.cdivnum i.year female pretrend shift posttrend [aw=pop] if year>2007 & year<=2016 & age!="All" & youth>0 & mcd=="Opioid", cluster(cdivnum)
lincom shift+posttrend*2
* (3f) youth interaction, then use full back period
reg deathrtscaled i.cdivnum i.year youth* female pretrend shift posttrend [aw=pop] if age!="All" & year>2007 & year<=2016 & mcd=="Opioid", cluster(cdivnum)
lincom shift+posttrend*2
lincom youthshift+youthpost*2
reg deathrtscaled i.cdivnum i.year youth* female pretrend shift posttrend [aw=pop] if age!="All" & year<=2016 & mcd=="Opioid", cluster(cdivnum)
outreg2 using "FigureA5A6A7.outreg.xls", `decs' `outstats' addtext("Substance","All opioids","Adults and minors","Separated","Years","1999-2016","Groups excluded","None") excel append
lincom shift+posttrend*2
lincom youthshift+youthpost*2

* (4) regression results: adding PDMP
reg deathrt i.cdivnum i.year female pretrend shift posttrend havePDMPmust [aw=pop] if year>2007 & year<=2013 & age=="All" & mcd=="T400T401T404", cluster(cdivnum)
lincom shift+posttrend*2
reg deathrt i.cdivnum i.year female pretrend shift posttrend havePDMPmust [aw=pop] if year<=2016 & age=="All" & mcd=="T400T401T404", cluster(cdivnum)
outreg2 using "FigureA5A6A7.outreg.xls", `decs' `outstats' addtext("Substance","Im opioids","Adults and minors","Pooled","Years","1999-2016","Groups excluded","None") excel append
lincom shift+posttrend*2
* (4b) youth interaction, switching mortality indicator to all opioid (from just Im), full back period
reg deathrtscaled i.cdivnum ib2010.year youth* PDMPmusty female pretrend shift posttrend havePDMPmust [aw=pop] if age!="All" & year<=2016 & mcd=="Opioid", cluster(cdivnum)
outreg2 using "FigureA5A6A7.outreg.xls", `decs' `outstats' addtext("Substance","All opioids","Adults and minors","Separated","Years","1999-2016","Groups excluded","None") excel append

* (5) predict "no regulation" hypotheticals
gen noregpoint=.
gen noregupper=.
gen noreglower=.
gen norefpoint=.
gen norefupper=.
gen noreflower=.
gen fitpoint=.
gen fitupper=.
gen fitlower=.
gen regimpactpoint=.
gen regimpactupper=.
gen regimpactlower=.
label var noregpoint "Adult - 25*child with no supply intervention"
label var norefpoint "Adult - 25*child with no reformulation effect"
label var fitpoint "Adult - 25*child fitted values"
* (5a) lincom
foreach yr of numlist 1999/2020 {
  lincom youth + youthoxy*usaoxy`yr' + youthpre*(`yr'-baseyear)*usaoxy`yr'
  replace noregpoint = -r(estimate) if year==`yr'                  // change child-adult to adult-child
  replace noregupper = noregpoint + 1.96*r(se) if year==`yr'
  replace noreglower = noregpoint - 1.96*r(se) if year==`yr'
  lincom youth + youthoxy*usaoxy`yr' + youthpre*(`yr'-baseyear)*usaoxy`yr' + PDMPmusty*usapdmp`yr'
  replace norefpoint = -r(estimate) if year==`yr'                  // change child-adult to adult-child
  replace norefupper = noregpoint + 1.96*r(se) if year==`yr'
  replace noreflower = noregpoint - 1.96*r(se) if year==`yr'
  lincom youth + youthoxy*usaoxy`yr' + youthpre*(`yr'-baseyear)*usaoxy`yr' + youthshift*0*usaoxy`yr' + youthpost*0*(`yr'-2011)*usaoxy`yr' + PDMPmusty*usapdmp`yr'
  replace fitpoint = -r(estimate) if year==`yr'                  // change child-adult to adult-child
  replace fitupper = fitpoint + 1.96*r(se) if year==`yr'
  replace fitlower = fitpoint - 1.96*r(se) if year==`yr'
}
foreach yr of numlist 1999/2010 {
  lincom havePDMPmust*usapdmp`yr' + PDMPmusty*usapdmp`yr'
  replace regimpactpoint = r(estimate) if year==`yr' & youth==1
  replace regimpactupper = regimpactpoint + 1.96*r(se) if year==`yr' & youth==1
  replace regimpactlower = regimpactpoint - 1.96*r(se) if year==`yr' & youth==1
  lincom havePDMPmust*usapdmp`yr'
  replace regimpactpoint = r(estimate) if year==`yr' & youth==0
  replace regimpactupper = regimpactpoint + 1.96*r(se) if year==`yr' & youth==0
  replace regimpactlower = regimpactpoint - 1.96*r(se) if year==`yr' & youth==0
}
foreach yr of numlist 2011/2020 {
  lincom youth + youthoxy*usaoxy`yr' + youthpre*(`yr'-baseyear)*usaoxy`yr' + youthshift*usaoxy`yr' + youthpost*(`yr'-2011)*usaoxy`yr' + PDMPmusty*usapdmp`yr'   // shift=1 because all of years are after 2010
  replace fitpoint = -r(estimate) if year==`yr'                  // change child-adult to adult-child
  replace fitupper = fitpoint + 1.96*r(se) if year==`yr'
  replace fitlower = fitpoint - 1.96*r(se) if year==`yr'
  lincom shift*usaoxy`yr' + posttrend*(`yr'-2011)*usaoxy`yr' + havePDMPmust*usapdmp`yr' + youthshift*usaoxy`yr' + youthpost*(`yr'-2011)*usaoxy`yr' + PDMPmusty*usapdmp`yr'   // shift=1 because all of years are after 2010
  replace regimpactpoint = r(estimate) if year==`yr' & youth==1
  replace regimpactupper = regimpactpoint + 1.96*r(se) if year==`yr' & youth==1
  replace regimpactlower = regimpactpoint - 1.96*r(se) if year==`yr' & youth==1
  lincom shift*usaoxy`yr' + posttrend*(`yr'-2011)*usaoxy`yr' + havePDMPmust*usapdmp`yr'   // shift=1 because all of years are after 2010
  replace regimpactpoint = r(estimate) if year==`yr' & youth==0
  replace regimpactupper = regimpactpoint + 1.96*r(se) if year==`yr' & youth==0
  replace regimpactlower = regimpactpoint - 1.96*r(se) if year==`yr' & youth==0
}
* (5b) chart
global cellfornationaldisplay `"mcd=="Opioid" & age=="18+" & female==0 & cdivnum==1 & year<=2016"'
#delimit ;
twoway (rcap noreglower noregupper year if $cellfornationaldisplay, lcolor(blue))
  (rcap fitlower fitupper year if $cellfornationaldisplay, lcolor(black))
  (scatter noregpoint year if $cellfornationaldisplay, msymbol(S) mcolor(blue) msize(small))
  (scatter fitpoint year if $cellfornationaldisplay, msymbol(S) mcolor(black) msize(small))
  (scatter navadultchildgap year if $cellfornationaldisplay, msymbol(O) mcolor(red) msize(small))
  (line noregpoint year if $cellfornationaldisplay, lcolor(blue))
  (line fitpoint year if $cellfornationaldisplay, lcolor(black)),
  title("Figure A5a.  Opioid fatality rates: adults vs. minors")
  subtitle("Adult rate - 25*(minor rate)")
  legend(order(4 "Fitted values" 3 "w/o supply intervention" 2 "95% confidence interval" 1 "95% confidence interval" 5 "Actual") region(lcolor(white)) size(small))
  graphregion(color(white)) ylabel(,nogrid) xtitle("") ytitle("Mortality gap")
  note("Note: Death rate for minors (ages 0-17) is multiplied by 25 to compare with adult rates (per 100K).  Squares and confidence" "intervals summarize estimates of the event-study model, zeroing out Post and PDMDP to simulate w/o supply intervention." "Population and mortality source: CDC Wonder.", size(vsmall))
  name(figurea5, replace);
#delimit cr
* (5c) impact chart by age
global cellfornationaldisplay `"mcd=="Opioid" & female==0 & cdivnum==1 & year>=2005 & year<=2016"'
#delimit ;
twoway (rcap regimpactlower regimpactupper year if $cellfornationaldisplay & age=="0-17", lcolor(blue))
  (rcap regimpactlower regimpactupper year if $cellfornationaldisplay & age=="18+", lcolor(black))
  (scatter regimpactpoint year if $cellfornationaldisplay & age=="0-17", msymbol(S) mcolor(blue) msize(small))
  (scatter regimpactpoint year if $cellfornationaldisplay & age=="18", msymbol(S) mcolor(black) msize(small))
  (line regimpactpoint year if $cellfornationaldisplay & age=="0-17", lcolor(blue))
  (line regimpactpoint year if $cellfornationaldisplay & age=="18+", lcolor(black)),
  title("Figure A5b.  Impact of the regulatory variables")
  subtitle("on opioid mortality rates, by age, with youth rates rescaled 25X")
  legend(order(4 "Ages 18+" 3 "Ages 0-17" 2 "95% confidence interval" 1 "95% confidence interval") region(lcolor(white)) size(small))
  graphregion(color(white)) ylabel(,nogrid) xtitle("") ytitle("Mortality per 100,000")
  note("Note: Death rate for minors (ages 0-17) is multiplied by 25 to compare with adult rates (per 100K).  Squares and confidence" "intervals summarize estimates of the combined effects of Post, Post-trend, and PDMP according to the event-study model." "Population and mortality source: CDC Wonder.", size(vsmall))
  name(figurea5b, replace);
#delimit cr
global cellfornationaldisplay `"mcd=="Opioid" & age=="18+" & female==0 & cdivnum==1 & year<=2016"'

* (6) Sensitivity analysis
* (6a) drop girls
reg deathrtscaled i.cdivnum ib2010.year youth* PDMPmusty female pretrend shift posttrend havePDMPmust [aw=pop] if age!="All" & (female==0 | age=="18+") & year<=2016 & mcd=="Opioid", cluster(cdivnum)
outreg2 using "FigureA5A6A7.outreg.xls", `decs' `outstats' addtext("Substance","All opioids","Adults and minors","Separated","Years","1999-2016","Groups excluded","Girls") excel append
* (6b) Log death rate, annual data (all but 3 zero-observation cells are girls)
reg lndeathrt i.cdivnum ib2010.year youth* PDMPmusty female pretrend shift posttrend havePDMPmust [aw=pop] if age!="All" & (female==0 | age=="18+") & year<=2016 & mcd=="Opioid", cluster(cdivnum)
outreg2 using "FigureA5A6A7.outreg.xls", `decs' `outstats' addtext("Substance","All opioids","Adults and minors","Separated","Years","1999-2016","Groups excluded","Girls") excel append
* (6c) Simulation for log death rate
*   (6c1) lincom
foreach yr of numlist 1999/2020 {
  lincom youth + youthoxy*usaoxy`yr' + youthpre*(`yr'-baseyear)*usaoxy`yr'
  replace noregpoint = -r(estimate) if year==`yr'                  // change child-adult to adult-child
  replace noregupper = noregpoint + 1.96*r(se) if year==`yr'
  replace noreglower = noregpoint - 1.96*r(se) if year==`yr'
  lincom youth + youthoxy*usaoxy`yr' + youthpre*(`yr'-baseyear)*usaoxy`yr' + PDMPmusty*usapdmp`yr'
  replace norefpoint = -r(estimate) if year==`yr'                  // change child-adult to adult-child
  replace norefupper = noregpoint + 1.96*r(se) if year==`yr'
  replace noreflower = noregpoint - 1.96*r(se) if year==`yr'
  lincom youth + youthoxy*usaoxy`yr' + youthpre*(`yr'-baseyear)*usaoxy`yr' + youthshift*0*usaoxy`yr' + youthpost*0*(`yr'-2011)*usaoxy`yr' + PDMPmusty*usapdmp`yr'
  replace fitpoint = -r(estimate) if year==`yr'                  // change child-adult to adult-child
  replace fitupper = fitpoint + 1.96*r(se) if year==`yr'
  replace fitlower = fitpoint - 1.96*r(se) if year==`yr'
}
foreach yr of numlist 2011/2020 {
  lincom youth + youthoxy*usaoxy`yr' + youthpre*(`yr'-baseyear)*usaoxy`yr' + youthshift*usaoxy`yr' + youthpost*(`yr'-2011)*usaoxy`yr' + PDMPmusty*usapdmp`yr'
  replace fitpoint = -r(estimate) if year==`yr'                  // change child-adult to adult-child
  replace fitupper = fitpoint + 1.96*r(se) if year==`yr'
  replace fitlower = fitpoint - 1.96*r(se) if year==`yr'
}
*   (6c2) chart
#delimit ;
twoway (rcap noreglower noregupper year if $cellfornationaldisplay, lcolor(blue))
  (rcap fitlower fitupper year if $cellfornationaldisplay, lcolor(black))
  (scatter noregpoint year if $cellfornationaldisplay, msymbol(S) mcolor(blue) msize(small))
  (scatter fitpoint year if $cellfornationaldisplay, msymbol(S) mcolor(black) msize(small))
  (scatter navlnadultchildgapnog year if $cellfornationaldisplay, msymbol(O) mcolor(red) msize(small))
  (line noregpoint year if $cellfornationaldisplay, lcolor(blue))
  (line fitpoint year if $cellfornationaldisplay, lcolor(black)),
  title("Figure A6.  Opioid fatality rates: adults vs. minors")
  subtitle("Log adult rate - Log minor rate")
  legend(order(4 "Fitted values" 3 "w/o supply intervention" 2 "95% confidence interval" 1 "95% confidence interval" 5 "Actual") region(lcolor(white)) size(small))
  graphregion(color(white)) ylabel(,nogrid) xtitle("") ytitle("Log mortality-rate gap")
  note("Note: Log specification, excluding girls.  All mortality rates are measured per 100K.  Squares and confidence" "intervals summarize estimates of the event-study model, zeroing out Post and PDMDP to simulate w/o supply intervention." "Population and mortality source: CDC Wonder.", size(vsmall))
  name(figurea6, replace);
#delimit cr

* (6d) Log death rate, biannual data (all but 3 zero-observation cells are girls)
*   (6d1) aggregate odd years with subsequent year
gen byte evenyr = (mod(year,2)==0)
replace year = year + 1 if evenyr==0
label var year "Current/lagged year if even/odd, respectively"
keep if age~="All" & mcd=="Opioid"
collapse (mean) deathrt youthpre youthpost youthshift youthoxy pretrend shift posttrend havePDMPmust PDMPmusty pop navlnadultchildgap, by(year cdivnum female youth)
*   (6d2) regression without girls
gen lndeathrt = log(deathrt)
reg lndeathrt i.cdivnum ib2010.year youth* PDMPmusty female pretrend shift posttrend havePDMPmust [aw=pop] if (female==0 | youth==0) & year<=2016, cluster(cdivnum)
outreg2 using "FigureA5A6A7.outreg.xls", `decs' `outstats' addtext("Substance","All opioids","Adults and minors","Separated","Years","1999-2016 biannual","Groups excluded","Girls") excel append
*   (6d2) regression with girls
reg lndeathrt i.cdivnum ib2010.year youth* PDMPmusty female pretrend shift posttrend havePDMPmust [aw=pop] if year<=2016, cluster(cdivnum)
outreg2 using "FigureA5A6A7.outreg.xls", `decs' `outstats' addtext("Substance","All opioids","Adults and minors","Separated","Years","1999-2016 biannual","Groups excluded","None") excel append
*   (6d3) lincom
gen noregpoint=.
gen noregupper=.
gen noreglower=.
gen norefpoint=.
gen norefupper=.
gen noreflower=.
gen fitpoint=.
gen fitupper=.
gen fitlower=.
forvalues yr = 2000(2)2020 {
  lincom youth + youthoxy*bianusaoxy`yr' + youthpre*(`yr'-0.5-baseyear)*bianusaoxy`yr'
  replace noregpoint = -r(estimate) if year==`yr'                  // change child-adult to adult-child
  replace noregupper = noregpoint + 1.96*r(se) if year==`yr'
  replace noreglower = noregpoint - 1.96*r(se) if year==`yr'
  lincom youth + youthoxy*bianusaoxy`yr' + youthpre*(`yr'-0.5-baseyear)*bianusaoxy`yr' + PDMPmusty*bianusapdmp`yr'
  replace norefpoint = -r(estimate) if year==`yr'                  // change child-adult to adult-child
  replace norefupper = noregpoint + 1.96*r(se) if year==`yr'
  replace noreflower = noregpoint - 1.96*r(se) if year==`yr'
  lincom youth + youthoxy*bianusaoxy`yr' + youthpre*(`yr'-0.5-baseyear)*bianusaoxy`yr' + youthshift*0*bianusaoxy`yr' + youthpost*0*(`yr'-0.5-2011)*bianusaoxy`yr' + PDMPmusty*bianusapdmp`yr'
  replace fitpoint = -r(estimate) if year==`yr'                  // change child-adult to adult-child
  replace fitupper = fitpoint + 1.96*r(se) if year==`yr'
  replace fitlower = fitpoint - 1.96*r(se) if year==`yr'
}
foreach yr of numlist 2011/2020 {
  lincom youth + youthoxy*usaoxy`yr' + youthpre*(`yr'-baseyear)*usaoxy`yr' + youthshift*usaoxy`yr' + youthpost*(`yr'-2011)*usaoxy`yr' + PDMPmusty*usapdmp`yr'
  replace fitpoint = -r(estimate) if year==`yr'                  // change child-adult to adult-child
  replace fitupper = fitpoint + 1.96*r(se) if year==`yr'
  replace fitlower = fitpoint - 1.96*r(se) if year==`yr'
}
*   (6d4) chart (girls are left out of the "actual" only because it is compared to fitted male cells
global bianndispcells `"female==0 & cdivnum==1 & year<=2016"'
#delimit ;
twoway (rcap noreglower noregupper year if $bianndispcells, lcolor(blue))
  (rcap fitlower fitupper year if $bianndispcells, lcolor(black))
  (scatter noregpoint year if $bianndispcells, msymbol(S) mcolor(blue) msize(small))
  (scatter fitpoint year if $bianndispcells, msymbol(S) mcolor(black) msize(small))
  (scatter navlnadultchildgapnog year if $bianndispcells, msymbol(O) mcolor(red) msize(small))
  (line noregpoint year if $bianndispcells, lcolor(blue))
  (line fitpoint year if $bianndispcells, lcolor(black)),
  title("Figure A7.  Opioid fatality rates: adults vs. minors")
  subtitle("Log adult rate - Log minor rate")
  legend(order(4 "Fitted values" 3 "w/o supply intervention" 2 "95% confidence interval" 1 "95% confidence interval" 5 "Actual") region(lcolor(white)) size(small))
  graphregion(color(white)) ylabel(,nogrid) xtitle("") ytitle("Log mortality-rate gap")
  note("Note: Biannual log specification, including girls.  All mortality rates are measured per 100K.  Squares and confidence" "intervals summarize estimates of the event-study model, zeroing out Post and PDMDP to simulate w/o supply intervention." "Population and mortality source: CDC Wonder.", size(vsmall))
  name(figurea7, replace);
#delimit cr



